# -*- coding: utf-8 -*-
"""
Created on Thu Oct 18 11:06:52 2012

@author: xudi
"""

# -*- coding: utf-8 -*-
"""
Created on Sun Oct 14 19:14:16 2012

@author: xudi
"""

import numpy as np
import pylab as pl
import copy
import gc

#physical constant define
c0=3e8;
eps0=8.85e-12
mu0=1.26e-6

#simulation constant
dz=1e-8
d_gra=1e-9
dt=dz/2/c0
sampling_rate=1/dt
Ntime=200000

f0=5e12
w0=2*np.pi*f0
t0=9*np.pi/(2*w0)
Tao=t0/4

#material constant uc0.13 
Npole=4
Pk=np.zeros(Npole,dtype=complex)
Pk[0]=-1.206064364027025
Pk[1]=2.00000000003219e12
Pk[2]=2.001995434351346e12-5.089190868348496e14j
Pk[3]=2.001995434351346e12+5.089190868348496e14j
#print Pk
Rk=np.zeros(Npole,dtype=complex)
Rk[0]=8.664242302240614e17*(d_gra/dz)
Rk[1]=-8.664020771160701e17*(d_gra/dz)
Rk[2]=(-6.116313254536545e12+1.559451294397289e15j)*(d_gra/dz)
Rk[3]=(-6.116313254536545e12-1.559451294397289e15j)*(d_gra/dz)
#print Rk


'''
#material constant uc0.13 
Npole=4
Pk=np.zeros(Npole,dtype=complex)
Pk[0]=-1.209216755621159
Pk[1]=2.00000000003219e12
Pk[2]=2.00199742308203e12+5.089190734129027e14j
Pk[3]=2.00199742308203e12-5.089190734129027e14j
#print Pk
Rk=np.zeros(Npole,dtype=complex)
Rk[0]=-8.664242302230176e17*(d_gra/dz)
Rk[1]=8.664020771142121e17*(d_gra/dz)
Rk[2]=(6.116229232747331e12-1.559451169231914e15j)*(d_gra/dz)
Rk[3]=(6.116229232747331e12+1.559451169231914e15j)*(d_gra/dz)
#print Rk
'''
'''
#material constant uc0.15
Npole=4
Pk=np.zeros(Npole,dtype=complex)
Pk[0]=-0.383079707748508
Pk[1]=2.00000000001016e12
Pk[2]=2.001486192995871e12+5.875137271947569e14j
Pk[3]=2.001486192995871e12-5.875137271947569e14j
#print Pk
Rk=np.zeros(Npole,dtype=complex)
Rk[0]=-9.98162323243687e17*(d_gra/dz)
Rk[1]=9.981431238417820e17*(d_gra/dz)
Rk[2]=(5.308934763269954e12-1.561859026197309e15j)*(d_gra/dz)
Rk[3]=(5.308934763269954e12+1.561859026197309e15j)*(d_gra/dz)
#print Rk
'''
'''
#material constant  uc0.2
Npole=4
Pk=np.zeros(Npole,dtype=complex)
Pk[0]=-0.035249513087656
Pk[1]=-2.000000000000103e12
Pk[2]=-2.000803000709099e12+7.838763956728979e14j
Pk[3]=-2.000803000709099e12-7.838763956728979e14j
#print 
Rk=np.zeros(Npole,dtype=complex)
Rk[0]=1.329639741650483e18*(d_gra/dz)
Rk[1]=-1.329625342058612e18*(d_gra/dz)
Rk[2]=(-3.98983129235818e12-1.56502124426222e15j)*(d_gra/dz)
Rk[3]=(-3.98983129235818e12+1.56502124426222e15j)*(d_gra/dz)
#print Rk
'''
'''
#material constant  uc0.4
Npole=4
Pk=np.zeros(Npole,dtype=complex)
Pk[0]=-0.002137281927179
Pk[1]=2.000000000000e12
Pk[2]=1.999054211752815e12+1.568790445808978e15j
Pk[3]=1.999054211752815e12-1.568790445808978e15j
#print Pk
Rk=np.zeros(Npole,dtype=complex)
Rk[0]=-2.658957852217624e18*(d_gra/dz)
Rk[1]=2.658950652402170e18*(d_gra/dz)
Rk[2]=(2.001070261580648e12-1.568144295330555e15j)*(d_gra/dz)
Rk[3]=(2.001070261580648e12+1.568144295330555e15j)*(d_gra/dz)
'''
alfak=(2-Pk*dt)/(2+Pk*dt)
print alfak
betak=2*Rk*eps0*dt/(2+Pk*dt)
print betak

#space size
nz=22000   #total grid numbers
size_SiO2=156
size_MgF=2000
sizefree=10000

#other infomation for graphene
#0.13
epsr_gra=5.948070532434306
#0.15
#epsr_gra=5.282858419236412
#0.2
#epsr_gra=4.206765308080218
#0.4
#epsr_gra=2.600731482016143
epsr_gramix=(dz-d_gra)/dz*eps0+(d_gra/dz)*epsr_gra
epsr_SiO2=3.9
epsr_MgF=5
gra_grids=1
gra_start=sizefree-gra_grids
SiO2_start=sizefree
SiO2_grids=size_SiO2
MgF_start=sizefree+size_SiO2
MgF_grids=size_MgF

#init memory 
Hy=np.zeros(nz)
Ex=np.zeros(nz)
R=np.zeros(Ntime)
T=np.zeros(Ntime)
#init the coefficient matrix
caEx=np.ones(nz)
cbEx=np.ones(nz)*dt/eps0
caHy=np.ones(nz)
cbHy=np.ones(nz)*dt/mu0

#modify the coefficient matrix in graphene region 
cbEx[SiO2_start:SiO2_start+SiO2_grids]=dt/eps0/epsr_SiO2
cbEx[MgF_start:MgF_start+MgF_grids]=dt/eps0/epsr_MgF
cbEx[gra_start:gra_start+gra_grids]=2*dt/(2*eps0*epsr_gramix+np.real(np.sum(betak)))
print cbEx[gra_start:gra_start+gra_grids]
print cbEx[SiO2_start:SiO2_start+SiO2_grids]
print cbEx[MgF_start:MgF_start+MgF_grids]
Jk=np.zeros((gra_grids,Npole),dtype=complex)
print "Good Luck"
print "test: the sum is the total current at each grid"
print np.sum(Jk,1)
#test a formula



#method used for update the field components
def updateHy():
    global Hy,Ex,caHy,cbHy,nz,dt
    Hy[1:nz-1]=caHy[1:nz-1]*Hy[1:nz-1]-cbHy[1:nz-1]*(Ex[1:nz-1]-Ex[:nz-2])/dz

def updateEx():
    global Ex,caEx,cbEx,nz,dt
    Ex[0:nz-1]=caEx[0:nz-1]*Ex[0:nz-1]-cbEx[0:nz-1]*((Hy[1:nz]-Hy[0:nz-1])/dz)
    
def updateJk(Ecurrent,Eold):
    global Jk,alfak,betak,dt,gra_grids
    for grid in range(0,gra_grids):
        Jk[grid]=alfak*Jk[grid]+betak*(Ecurrent[grid]-Eold[grid])/dt

def updateGrapheneEx(Zstart,width):
    global Ex, cbEx
    sumJk=np.real(np.sum((alfak+1)*Jk/2.0,1))
    Ex[Zstart:Zstart+width]=Ex[Zstart:Zstart+width]+cbEx[Zstart:Zstart+width]*sumJk
    
def source(realtime):
    G=-1*np.cos(w0*realtime)*np.exp(-(realtime-t0)**2/Tao**2)
    return G

for qtime in range(0, Ntime):
    UHy0=Hy[0]
    UHy1=Hy[1]
    UHy2=Hy[nz-1]
    UHy3=Hy[nz-2]
    
    updateHy()
    
    Hy[0]=UHy1+((c0*dt-dz)/(c0*dt+dz))*(Hy[1]-UHy0)
    Hy[nz-1]=UHy3+((c0*dt-dz)/(c0*dt+dz))*(Hy[nz-2]-UHy2)
      
    Hy[5000]+=(dt/(mu0*dz))*source(qtime*dt)
    
    Eold=copy.copy(Ex[gra_start:gra_start+gra_grids])
    
    updateEx()
    
    updateGrapheneEx(gra_start,gra_grids)
    
    Ecurrent=copy.copy(Ex[gra_start:gra_start+gra_grids])
    updateJk(Ecurrent,Eold)
    
    Ex[5000]+=((dt/(eps0*dz))*(eps0/mu0)**0.5)*source((qtime+1.5)*dt)
    
    R[qtime]=Ex[4500]
    T[qtime]=Ex[21000]

#np.savetxt('Reflection(time)',R)

#plot source R T and RFFT TFFT
t=np.arange(0,Ntime*dt,dt)
Xs=source(t)

fft1=np.fft.fft(T,10000000)
fft2=np.abs(np.fft.fftshift(fft1))
fft3=np.fft.fft(R,10000000)
fft4=np.abs(np.fft.fftshift(fft3))
fft5=np.fft.fft(Xs,10000000)
fft6=np.abs(np.fft.fftshift(fft5))
freqs=np.linspace(0,sampling_rate,num=10000000,endpoint=False)-sampling_rate/2
y1=np.arange(0,Ntime*dt,dt)


pl.figure(figsize=(8,9))

pl.subplot(311)
pl.plot(y1,T,'r')
pl.plot(y1,R,'b')
pl.ylim(-1,1)
pl.xlabel("Time (s)")
pl.ylabel("Gauss inedx")
pl.legend((u"Tr ","Rr "))

pl.subplot(312)
pl.xlim(1e12,1e13)
pl.plot(freqs,fft2,'r')
pl.plot(freqs,fft4,'b')
pl.xlabel(" frequence (Hz)")
pl.ylabel("Reflection and Transmission inedx")
pl.legend((u"Tr (fft)","Rr (fft)"))

dbvalue=20*np.log10(fft4/fft6)
#np.save("Muc0.22",dbvalue)

pl.subplot(313)
pl.xlim(2e12,4e12)
pl.ylim(-50,0)
#pl.plot(freqs,fft2/fft6,'r')
pl.plot(freqs,dbvalue,'b')
pl.xlabel(" frequence (Hz)")
pl.ylabel("Refraction and Transmission scale inedx")
pl.legend((u"Rr/Gauss "))
pl.grid('on')  



